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ABSTRACT 


The  Source  Physics  Experiment  (SPE-N)  at  the  Nevada  National  Security  Site  is  planned  as  a  series  of  chemical 
explosions  under  a  variety  of  emplacement  conditions.  The  goal  of  the  SPE-N  is  to  improve  our  physical 
understanding  and  ability  to  model  explosively  generated  seismic  waves,  particularly  S-waves.  The  first  SPE 
explosion  (SPE1)  consisted  of  a  100  kg  shot  at  a  depth  of  60  meters  in  granite  (Climax  Stock).  The  shot  was  well- 
recorded  by  an  array  of  over  150  instruments,  including  both  near- field  wave  motion  measurements  as  well  as  far- 
field  seismic  measurements.  This  paper  focuses  on  measurements  and  modeling  of  the  near- field  data.  A 
complimentary  study  focusing  on  far- field  seismic  observations  is  the  subject  of  a  companion  paper  (Mellors  et  al., 
these  Proceedings).  The  near- field  data  included  triaxial  acceleration  measurements  at  eighteen  different  locations 
azimuthally  distributed  around  the  explosive  charge.  Three  triaxial  accelerometers  were  embedded  in  each  of  six 
vertical  boreholes,  distributed  in  two  concentric  rings  around  the  charge.  The  inner  ring  consisted  of  three 
equidistant  boreholes  at  a  radius  of  1 0  m  from  the  charge,  and  the  outer  ring  consisted  of  another  three  equidistant 
boreholes  at  a  radius  of  20  m.  In  each  borehole,  the  accelerometers  were  vertically  distributed  at  depths  of  60  m 
(shot  horizon),  50  m  and  15  m.  Surface  accelerations  were  also  recorded  along  a  radial  line  centered  at  surface 
ground  zero.  A  review  of  the  SPE1  data  shows  that  the  peak  radial  velocity  as  a  function  of  scaled  range  is 
consistent  with  previous  nuclear  explosion  data  but  exhibits  greater  variability.  The  scaled  peak  radial  displacement 
variation  also  exhibits  greater  variability  but  the  mean  values  are  significantly  higher  than  exhibited  in  previous 
nuclear  explosion  data.  These  higher  displacements  were  also  observed  in  calculations  performed  with  an  ensemble 
constitutive  model,  which  was  based  on  nuclear  explosion  data  in  hard  rock,  but  employed  a  Jones-Wilkins-Lee 
(JWL)  equation  of  state  for  the  ammonium  nitrate  fuel  oil  (ANFO)  explosive  mix  used  in  SPE1.  The  reason  for  this 
behavior  is  believed  to  be  the  higher  effective  ratio  of  specific  heats  in  the  explosion  products  of  the  chemical 
explosive,  leading  to  higher  residual  cavity  pressure.  Azimuthal  scatter  in  the  velocity  data  correspond  to  joint 
orientation,  and  as  anticipated,  the  joints  appear  to  be  the  principal  source  of  observed  shear  wave  generation  in  the 
near- field.  Preliminary  modeling  of  the  SPE1  data  shows  that  continuum  simulations  that  do  not  explicitly  account 
for  the  effect  of  joints  will  not  successfully  reproduce  the  observed  directional  variations  in  the  recorded  data. 
However,  2D  and  3D  simulations  that  explicitly  account  for  joints  and  pre-existing  fractures  show  that  a  low  friction 
angle,  derived  with  water-filled  joints,  may  account  for  the  observed  variation  in  peak  velocity  and  displacement. 
Waves  appear  to  propagate  more  readily  in  the  direction  of  persistent  joints,  as  opposed  to  staggered  joints. 
Furthermore,  the  anisotropy  associated  with  wave  propagation  seems  to  be  more  pronounced  when  the  friction  angle 
was  lowered  to  account  for  the  effect  of  saturation.  Further  modeling  will  be  conducted  with  continued  focus  on  the 
effect  of  the  presence  of  joints  and  their  properties  on  shear  wave  generation.  Simulation  results  will  be  compared  to 
experimental  measurements  of  both  radial  and  non-radial  motions  from  the  SPE1  event,  as  well  as  from  planned 
future  SPE  explosions. 


Note:  The  Source  Physics  Experiments  at  the  Nevada  Nuclear  Security  Site  (SPE-N)  in  2011  should  not  be  confused 
with  the  2003  Source  Phenomenology  Experiments  conducted  in  Arizona  (SPE-A)  (Yang  and  Bonner,  2009). 
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OBJECTIVES 


The  SPE-N  is  planned  as  a  series  of  chemical  explosions  under  a  variety  of  emplacement  conditions.  The  SPE-N 
goal  is  to  improve  our  physical  understanding  and  ability  to  model  how  explosions  generate  seismic  waves, 
particularly  S-waves. 

RESEARCH  ACCOMPLISHED 


Continuum  Modeling 

An  Eulerian  hydrocode,  GEODYN,  was  used  to  calculate  spherical  wave  propagation  both  in  2D  and  3D.  The 
(ensemble)  constitutive  model  employed  was  calibrated  to  match  peak  velocities  and  displacements  for  a  number  of 
underground  nuclear  explosions  of  different  magnitude  in  granitic  rock  conducted  both  in  the  US  and  abroad 
(Antoun  et  al.,  2002). 


Figure  1.  Peak  velocity  and  displacement  calculated  with  the  ensemble  model  compared  with  SPE1  data.  The 
black  lines  are  least  squares  power-law  fits  to  the  ensemble  constitutive  model  with  nuclear  source 
and  the  red  lines  are  least  squares  power-law  fits  to  the  ensemble  model  with  SPE-N  ANFO  source. 

Figure  1  shows  that  there  is  only  a  slight  difference  between  the  calculated  velocity  attenuation  with  nuclear  or 
chemical  sources  and  both  represent  the  SPE1  data  reasonably  well;  though  the  latter  exhibit  considerable  scatter. 

On  the  other  hand,  the  calculated  displacement  is  higher  by  roughly  a  factor  of  2  for  an  ANFO  explosion  compared 
to  a  nuclear  explosion  of  the  same  yield.  This  trend  is  also  exhibited  by  the  measured  data.  The  reason  for  this 
behavior  is  believed  to  be  the  higher  effective  ratio  of  specific  heats  in  the  explosion  products  of  the  chemical 
explosive,  leading  to  higher  residual  cavity  pressure  (Glenn  and  Goldstein,  1994). 

The  ensemble  model  is  an  isotropic  plasticity  model  designed  for  large  scale  simulations.  It  does  not  use  small  scale 
test  data  for  calibration  and  assumes  that  any  joints  in  the  granite  formation  are  distributed  uniformly  resulting  in  an 
isotropic  material  response.  Since  the  model  is  isotropic,  it  cannot  explain  the  observed  variability  in  the  data, 
deviations  from  spherical  wave  propagation  or  the  generation  of  significant  shear  (transverse)  waves  from 
explosions.  However,  even  the  closest  gage  location  at  10  m  indicates  the  wave  motion  data  from  SPE1  exhibits 
significant  non-radial  motions.  These  non-radial  motions  are  not  symmetrically  distributed  around  the  source,  and 
like  the  radial  motions,  they  exhibit  significant  azimuthal  variations.  Data  from  the  instrumentation  bore  holes,  as 
well  as  from  the  SPE1  charge  emplacement  hole  indicate  that  the  granite  is  highly  fractured,  with  an  average  joint 
spacing  of  about  1ft  (30.5  cm).  In  addition,  the  data  shows  that  the  joints  are  not  randomly  oriented;  instead,  four 
joint  sets  have  been  identified;  each  of  which  exhibit  a  distinct  orientation  with  relatively  little  variability.  The  high 
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angle  joint  sets  identified  from  the  data  are  shown  in  Figure  2(a).  Figure  2(b)  is  a  graphic  depiction  of  the  peak 
radial  velocities  recorded  at  all  gage  locations  in  the  SPE1  near- field.  The  vectors  shown  are  oriented  in  the  direction 
of  the  sensor  and  their  magnitudes  are  proportional  to  the  magnitude  of  the  measured  peak  radial  velocity.  For  ease 
of  presentation,  they  have  all  been  aggregated  around  the  circumference  of  a  circle  centered  on  the  charge.  This 
figure  appears  to  show  a  correlation  between  the  orientation  of  the  joints  and  the  preferential  directions  of  energy 
radiation.  Even  though  this  correlation  is  not  conclusive  due  to  the  loss  of  data  at  some  of  the  gage  stations,  it  is 
sufficiently  compelling  to  advance  the  hypothesis  that  the  rock  fabric,  as  manifested  by  the  joints  and  fractures,  is 
primarily  responsible  for  the  observed  non-radial  motions  as  well  as  the  azimuthal  variations  in  the  near-field 
motion  measurements.  The  remainder  of  this  paper  will  examine  this  hypothesis  utilizing  detailed  mesoscale 
simulations  that  discretely  incorporate  the  joints  and  account  for  slip  along  their  surfaces. 


(a)  Dip  direction  for  high  anglefractures  (b)  Scaled  peak  radial  velocity  measurements 

Figure  2.  Correlation  between  the  orientation  of  the  joint  sets  in  the  SPE1  test  bed  (a)  and  the  radial  motions 
recorded  at  all  the  near-field  sensor  locations. 

Discrete  Modeling 

In  the  discrete  mesoscale  simulations,  the  geologic  medium  is  represented  as  an  assembly  of  intact  rock  blocks 
separated  by  joints  and  fractures.  The  intact  rock  is  modeled  using  a  constitutive  model  that  accounts  for  the  salient 
features  of  the  response  of  geologic  materials,  and  calibrated  for  granite  using  laboratory  data,  including  data  from 
core  samples  obtained  from  the  SPE1  test  bed.  The  joints  are  modeled  as  frictional  interfaces  with  their  own 
constitutive  relations,  which  include  governing  equations  for  both  frictional  sliding  along  the  joint  surface  as  well  as 
compressibility,  or  compliance,  normal  to  the  joint  surface.  The  computational  resources  required  in  this  case  are 
formidable  because  of  the  large  number  of  joints  and  the  3D  character  of  the  problem.  Preliminary  analysis  done  in 
2D  showed  that  the  joint  orientation  and  persistence  play  a  big  role  in  formation  and  attenuation  of  the  waves  (see 
Figure  3). 
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Figure  3.  Pressure  contours  (0-200  MPa)  calculated  in  a  2-set  joint  system  with  different  persistency  (A- 
100%,  B-75%  and  C-25%). 

By  decreasing  the  joint  persistency  from  100%  (Figure  3(a))  to  below  75%  (Figure  3(b)  and  (c))  an  anisotropic  wave 
becomes  a  symmetric  wave  that  would  be  generated  in  an  isotropic  media. 

For  the  SPE1  calibration  shot,  we  used  a  Lagrangian  hydrocode,  GEODYN-L  (Vorobiev,  2010)  to  model  3D  wave 
propagation  in  a  region  100x100x100  m  with  two  joint  spacings,  3  m  and  1.5  m  to  study  the  effect  of  joints  on  the 
wave  propagation. 


Granite  Model  for  Small  Samples  and  its  Extension  to  Large  Scale 


The  model  developed  in  (Liu  et.  al.,  2005;  Vorobiev  et.  al.,  2007;  Vorobiev,  2008)  was  used  to  describe  the 
constitutive  behavior  of  the  intact  granite.  It  describes  the  effects  of  bulking,  porous  compaction,  shear  enhanced 
compaction,  and  pressure  dependent  yield  on  the  material  response.  The  yield  strength  function  is  chosen  in  the  form 
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where  (3  is  the  Lode  angle,  Yc  is  the  unconfined  compressive  strength,  8  is  a  function  of  plastic  strain  describing 
material  hardening,  Pc  and  P0  are  history  dependent  variables  describing  a  cap  surface  that  governs  porous 

compaction  and  Ym  is  the  ultimate  compressive  strength.  The  coefficients  s  and  m  depend  on  the  quality  of  the 

rock  as  expressed  by  the  Geological  Strength  Index  (GSI),  a  rock  quality  index  widely  used  by  the  geotechnical 
engineering  community: 
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For  the  intact  material,  s=  1  and  the  value  mi  can  be  found  from  static  laboratory  tests  on  intact  samples. 


(2) 


417 


2011  Monitoring  Research  Review:  Ground-Based  Nuclear  Explosion  Monitoring  Technologies 


Figure  4.  Comparison  of  model  fits  to  unconfined  compressive  strength  test  for  intact  granite  samples  from 
the  SPE-N  test  bed.  Bold  lines  are  simulation  results,  and  thinner  lines  with  markers  are 
experimental  data  (Broome,  2011). 

The  model  described  by  Equations  (1)  and  (2)  was  calibrated  to  match  intact  granite  data  from  the  literature,  as  well 
as  the  unconfmed  compressive  tests  performed  on  the  granite  samples  collected  at  the  SPE1  site  (see  Figure  4).  The 
model  can  be  scaled  to  field  scale  by  changing  the  GSI  index,  introducing  extra  porosity  associated  with  the  joints 
and  modifying  the  porous  compaction  law.  Figure  5  shows,  that  while  the  small  scale  model  does  not  reproduce  the 
wave  attenuation  data  for  large  scale  explosions,  such  as  SPE1,  it  can  be  tuned  to  match  the  data  by  changing  the 
GSI  index,  which  describes  the  quality  of  rock  mass.  In  this  case,  a  GSI=60  was  needed  to  match  the  measured 
displacement,  which  appropriately  corresponds  to  a  blocky  rock  with  fair-to-poor  surface  conditions.  In  the 
mesoscale  simulations  described  below,  we  used  the  model  for  intact  granite,  because  the  joints  are  included 
explicitly  in  the  simulations. 

Model  for  Dry  and  Wet  Joints  Informed  by  Meso-scale  Modeling 

One  complicating  factor  in  performing  the  discrete  mesoscale  simulations  described  here  is  the  lack  of  joint 
characterization  data  to  calibrate  the  contact  parameters  used  in  the  simulations.  This  is  further  exacerbated  by  the 
presence  of  a  perched  water  table  at  the  location  of  the  SPE1  shot,  thus  introducing  added  complications  to  the 
behavior  of  the  joints.  To  alleviate  this  shortcoming,  and  provide  the  information  needed  for  a  reasonable 
representation  of  the  response  of  the  joints  in  the  mesoscale  simulations,  a  series  of  detailed  simulations  were 
performed  to  examine  the  behavior  of  a  single  joint  under  both  dry  and  saturated  conditions.  Hydrocode  simulations 
are  performed  here  with  GEODYN,  which  is  a  massively  parallel  Eulerian  hydrocode  with  AMR  capabilities.  The 
“synthetic  data”  obtained  from  the  single  joint  simulations  were  used  to  calibrate  the  joint  model  utilized  in  the 
mesoscale  simulations  of  the  SPE1  near- field. 
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Figure  5.  Peak  displacement  and  velocity  attenuations  calculated  in  ID  spherical  explosion  using  the  small 

scale  model  (solid  squares  and  circles)  versus  the  nuclear  data  fit.  Open  squares  and  circles  show  an 
in  situ  model  calculation. 


The  maximum  shear  strength  of  a  rock  joint  depends  on  the  normal  stress.  Coulomb’s  linear  relationship  describes 
this  dependence  well  within  the  range  of  experimental  data: 

T  =  c  +  an  tan^  (3) 


where  ris  the  peak  shear  stress  crn  is  the  normal  stress,  c  is  the  cohesion  intercept  and  ^  is  the  friction  angle.  It  is 
assumed  that  these  stresses  are  averaged  along  the  length  of  the  crack.  A  wider  ranging  empirical  model  was 
introduced  by  Barton  to  predict  the  maximum  possible  shear  for  a  given  applied  normal  stress  that  accounts  for 
curvature  in  the  strength  envelopes  for  non-planar  rock  joints  (Barton,  1973;  Barton  and  Choubey,  1977).  Barton 
introduces  the  joint  roughness  coefficient  (, JRC ),  the  joint  compressive  strength  (. JCS)  and  a  basic  friction  angle  ($,) 
in  his  model,  which  have  made  it  quite  popular  due  to  the  simple  determination  of  these  parameters.  For  example, 
the  JRC  is  determined  by  visual  inspection.  In  this  section  we  characterize  the  roughness  of  a  simulated  joint  using 
the  ideas  from  Barton  and  compare  a  dry  and  wet  joint  for  the  granite  at  the  SPE1  site. 


To  characterize  the  joint  we  simulate  a  joint  surface  with  a  specific  roughness  that  is  correlated  to  the  commonly 
used  JRC  values.  The  successive  random  midpoint  method  is  used  to  perturb  a  500  mm  line  as  shown  in  Figure  6. 
To  objectively  quantify  the  JRC  value,  the  fractal  dimension  must  be  related  through  an  empirical  expression  (Lee 
et.  al.,  1990): 


JRC  =  -0.87804  +  37.7844 


"  D  -  I 
V0.015; 


16.9304 


f  Dj- lj2 

'vO.Ol  5y 


(4) 


where  D  is  the  fractal  dimension.  The  fractal  dimension  of  the  line  shown  in  Figure  6  is  D=  1.00548  which 
corresponds  to  JRC  =  11  from  Equation  (4). 
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Figure  6.  GEODYN  simulation  and  results,  (a)  Part  of  the  500  mm  long  joint  profile  is  shown  at  left  with 

damage  indicated  by  the  color  overlay.  The  left  side  of  the  fractal  joint  is  moving  upward  at  5  m/s 
and  the  right  is  moving  downward  with  the  same  speed.  At  the  end  of  the  simulation  the  stresses  in 
each  vertical  slice  from  the  boxed  region  are  integrated  for  every  time  step,  (b)  The  integrated 
normal  and  shear  stress  pairs  with  a  least  squares  fit  to  Equation  (3). 

Figure  6  shows  the  results  for  8  different  simulations  for  a  dry  and  wet  joint  with  a  mid-range  roughness  coefficient 
(JRC=1 1).  The  left  side  of  the  fractal  joint  is  moving  upward  at  5  m/s  and  the  right  is  moving  downward  with  the 
same  speed.  Initially,  a  linear  vertical  velocity  gradient  is  specified  across  the  joint.  The  initial  gap  between  the  two 
matched  surfaces  is  2.5  mm.  At  the  end  of  the  simulation  the  stresses  in  each  vertical  slice  from  the  boxed  region  are 
integrated  for  every  time  step.  The  integrated  normal  and  shear  stress  pairs  are  shown  in  Figure  6(b)  and  the 
maximum  shear  stress  is  determined  and  plotted  in  the  figure  at  right.  A  least  squares  fit  of  this  data  produces  the 
cohesion  and  averaged  friction  coefficient  (see  Equation  3)  for  the  velocity  range  detected  at  10  m  distance  in  the 
SPE1  experiment.  The  values  for  the  black  and  red  line  in  the  figure  corresponding  to  the  two  terms  in  Equation  (3) 
are:  (dry)  c  =1.74,  tan($=0.67,  (wet)  c  =1.24,  tan(^)=0.21. 

Discrete  Problem  Setup 

The  jointed  rock  formation  was  created  by  paving  the  space  with  parallelepiped  blocks.  Each  block  was 
subdiscretized  by  a  number  of  elements  with  contacts  set  at  the  exterior  faces.  The  properties  of  the  contacts  were 
varied  randomly  from  sticky  to  weak  contact  (with  low  frictional  coefficient,  corresponding  to  the  wet  joint 
described  above).  In  addition  to  the  frictional  properties,  the  persistence  of  the  joints  was  controlled  by  staggering 
the  blocks  in  two  directions.  This  approach  allowed  us  to  describe  a  joint  system  with  three  main  joint  directions. 
These  directions  were  identified  from  the  experimental  characterization  of  the  joints  at  different  locations.  Figure  7 
displays  the  3D  problem  setup  and  Figure  8  shows  the  measured  joint  dip  directions,  as  well  as  the  orientation  of  the 
XY  plane,  used  in  the  calculations.  The  X  axis  was  oriented  with  joint  set  #1  (see  Figure  2(a)  for  joint  set 
definition).  Note  that  we  have  treated  the  two  closely  oriented  joint  sets  #1  and  #2  from  Figure  2(a)  and  Figure  8  as  a 
one  joint  set  #1  in  the  simulations  as  indicated  by  the  red  arrow  in  Figure  7  (a). 
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Figure  7.  3D  problem  initialization  for  rock  formation  with  3  joint  sets.  The  Z  axis  is  oriented  towards  the 
surface.  The  blocks  around  the  source  (at  depth  Z=-55  m)  were  further  refined. 


R=10 


Figure  8.  Joint  dip  directions,  well  locations  and  coordinate  axis  locations. 


Simulation  results 

The  results  of  our  discrete  simulations  indicate  that  joints  may  induce  a  large  amount  of  anisotropy.  Figure  9  shows 
that  there  is  a  significant  variation  in  attenuation  of  radial  velocity  in  the  horizontal  plane  at  1 80  ft  (explosion  depth) 
at  ranges  of  10  and  20  m  and  at  ranges  of  14  m  and  22  m  at  a  depth  of  150  ft  (10  m  above).  This  can  be  explained  by 
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sliding  that  occurs  at  the  joints  and  is  confirmed  by  the  discrete  calculations.  Both  the  discrete  calculations 
performed  with  GEODYN-L  and  the  experiments  show  significant  scatter  in  peak  velocity  and  peak  displacement 
values  measured  at  the  same  range  at  different  azimuths,  which  indicates  that  the  joints  have  a  profound  effect  on 
the  wave  anisotropy  and  can  help  explain  the  experimental  observations.  Figure  10  shows  radial  velocity  histories 
calculated  at  two  ranges,  10  m  and  20  m.  As  in  the  experiment,  the  peak  velocities  are  different  at  various  azimuths. 
The  velocities  shown  in  Figure  10  correspond  to  X  and  Y  directions,  which  are  the  most  favorable  for  wave 
propagation.  Similar  to  the  calculations,  the  highest  velocity  was  measured  at  location  #1,  which  is  shown  with  the 
solid  line.  Results  of  the  simulations  underestimate  both  the  velocity  and  displacement  at  location  #1.  This  is  likely 
due  to  the  fact  that  we  have  used  roughly  3  times  larger  joint  spacing  than  extant  from  observations.  The  calculation 
cost  increases  rapidly  with  the  number  of  joints  treated  by  the  contacts  and  there  were  insufficient  resources 
available  in  our  preliminary  study.  Code  improvements  are  underway  to  enable  more  realistic  simulations  that  utilize 
the  actual  joint  spacing.  Preliminary  results  indicate  that  both  peak  velocity  and  the  displacement  should  increase 
when  more  joints  are  added.  Moreover,  the  results  show  that  taking  the  joints  into  consideration  dramatically 
changes  the  wave  forms. 

The  calculations  shown  in  Figure  10  were  done  with  the  intact  model  without  joints  (GEODYN),  as  well  as  with  the 
jointed  model  (GEODYN-F).  The  GEODYN  calculations  produced  peak  velocities  a  factor  of  2,  and  peak 
displacements  ~20  times,  smaller  than  those  observed  in  the  SPE1  experiment,  whereas  the  GEODYN-E  simulations 
that  included  the  joints  produced  results  that  are  comparable  to  the  data,  once  again  highlighting  the  important  role 
that  the  rock  fabric  plays  in  shaping  the  waves  produced  by  subsurface  explosions. 


Figure  9.  Peak  velocity  and  displacement  versus  range  measured  in  experiment  (crosses)  and  calculated  using 
discrete  model  (circles  and  squares).  The  solid  line  shows  the  nuclear  explosion  data  fit. 
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Figure  10.  Radial  velocity  evolution  measured  in  experiments  and  calculated  with  discrete  modeling  (green 

and  red  curves)  using  1.5  m  joint  spacing  and  12  degree  friction  angle  at  two  different  ranges,  10  m 
and  20  m.  Comparison  is  also  made  with  continuum  calculations  with  the  intact  model  without 
joints. 


CONCLUSIONS  AND  RECOMMENDATIONS 


Results  from  SPE1  show  peak  radial  velocity  as  a  function  of  scaled  range  is  consistent  with  previous  nuclear 
explosion  data  but  exhibits  somewhat  greater  variability.  The  scaled  peak  radial  displacement  variation  also  exhibits 
greater  variability  but  the  mean  values  are  significantly  higher  than  exhibited  in  previous  nuclear  explosion  data. 

The  reason  for  the  larger  displacements  is  believed  to  be  the  higher  effective  ratio  of  specific  heats  in  the  explosion 
products  of  the  chemical  explosive,  leading  to  higher  residual  cavity  pressure.  Preliminary  2D  and  3D  calculations 
suggest  that  the  observed  variations  in  both  peak  velocity  and  displacement  can  be  explained  by  the  presence  of 
water-filled  joints.  We  have  shown  that  using  discrete  simulations,  where  joints  in  the  rock  formation  are  treated 
explicitly,  helps  to  understand  the  anisotropy  in  wave  propagation.  The  simulated  velocity  and  displacement 
histories  were  sensitive  to  average  joint  spacing,  statistical  variation  in  joint  spacing,  and  joint  persistency.  Site 
specific  information  on  these  joint  properties  would  better  constrain  the  models  and  improve  the  quality  of  the 
results. 

Finally,  it  should  be  noted  that  in  the  discrete  simulations  we  used  a  constitutive  model  for  the  rock  that  was 
obtained  by  calibration  with  laboratory-scale  rock  samples  collected  from  the  site  which  is  in  contrast  to  traditional 
approaches  where  laboratory-scale  data  are  arbitrarily  modified  to  fit  field-scale  results. 
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